home *** CD-ROM | disk | FTP | other *** search
/ PsL Monthly 1993 December / PSL Monthly Shareware CD-ROM (December 1993).iso / prgmming / dos / pascal / p_mat.exe / TESTOP.PAS < prev   
Encoding:
Pascal/Delphi Source File  |  1993-01-17  |  7.7 KB  |  280 lines
  1. Program testop;
  2. { test suite for pmatop }
  3.  
  4. Uses pmat, pmatop;
  5.  
  6. Procedure testMfunctions;
  7. Var
  8.   x,y : vmatrixptr;
  9. Begin
  10.     Dispatch^.Inclevel;
  11.     new( x, makematrix( 1, 1 ) );
  12.     new( y, makematrix( 1, 1 ) );
  13.     
  14.     x := matequals( x, add( Ident( 5 ), fill( 5, 5, 1 ) ) );
  15.     { sort rows in ascending order of the first column }
  16.     y := matequals( y, MSort( x, 1, 1 ) );
  17.     y^.Show( 'sorted on rows' );
  18.     
  19.     { sort columns in ascending order of the second row }
  20.     y := matequals( y, MSort( y, 2, 0 ) );
  21.     y^.Show( 'sorted on columns' );
  22.     
  23.     y := matequals( y, Mexp( x ) );
  24.     y^.Show( 'exponent of elements' );
  25.     
  26.     y := matequals( y, Mabs( x ) );
  27.     y^.Show( 'absolute value of elements' );
  28.     
  29.     y := matequals( y, Mlog( x ) );
  30.     y^.Show( 'natural logrithm of elements' );
  31.     
  32.     y := matequals( y, Msqrt( x ) );
  33.     y^.Show( 'square root of elements' );
  34.     
  35.     dispose( x, killvmatrix );
  36.     dispose( y, killvmatrix );
  37.     Dispatch^.Declevel;
  38. End;
  39.  
  40. Procedure testPatterned;
  41. Var
  42.   H, K : vmatrixptr;
  43. Begin
  44.     Dispatch^.Inclevel;
  45.     new( H, makematrix( 5, 5 ) );
  46.     new( K, makematrix( 5, 5 ) );
  47.     
  48.     h := matequals( h, helm( 5 ) );
  49.     h^.Show( 'Helmert Matrix' );
  50.     writeln;
  51.     writeln( 'trace of the helmert(5): ', trace( H ): 7: 5 );
  52.     writeln( 'Determinant of Helmert Matrix: ', Det( H ): 7: 5 );
  53.     
  54.     k := matequals( k, index( 1, 5, 1 ) );
  55.     k^.Show( '1 5 1' );
  56.     k := matequals( k, index( 5, 1, - 1 ) );
  57.     k^.Show( '5 1 -1' );
  58.     
  59.     k := matequals( k, index( 1, 5, 2 ) );
  60.     k^.Show( '1 5 2' );
  61.     k := matequals( k, index( 5, 1, - 2 ) );
  62.     k^.Show( '5 1 -2' );
  63.     
  64.     k := matequals( k, Kron( Ident( 2 ), H ) );
  65.     k^.Show( 'I@H' );
  66.     
  67.     k := matequals( k, H );
  68.     k^.mm( 1, 1 )^ := 1.0E - 9;
  69.     writeln;
  70.     writeln( 'Determinant of Pivoted Matrix: ', Det( k ): 7: 5 );
  71.     
  72.     dispose( h, killvmatrix );
  73.     dispose( k, killvmatrix );
  74.     Dispatch^.Declevel;
  75. End;
  76.  
  77. Procedure testDecomp;
  78. Var
  79.   H, S, Q, D, R : vmatrixptr;
  80. Begin
  81.     Dispatch^.Inclevel;
  82.     new( H, makematrix( 5, 5 ) );
  83.     new( S, makematrix( 5, 5 ) );
  84.     new( Q, makematrix( 5, 5 ) );
  85.     new( D, makematrix( 5, 5 ) );
  86.     new( R, makematrix( 5, 5 ) );
  87.     
  88.     h := matequals( h, Add( Ident( 5 ), Fill( 5, 5, 1 ) ) );
  89.     S := matequals( S, Cholesky( h ) );
  90.     S^.Show( 'Cholesky decomposition' );
  91.     S := matequals( S, mult( Tran( S ), S ) );
  92.     S^.Show( 'S`S' );
  93.     
  94.     QR( h, Q, R, true );
  95.     Q^.Show( 'Q' );
  96.     R^.Show( 'R' );
  97.     S := matequals( S, mult( Q, R ) );
  98.     S^.Show( 'QR' );
  99.     
  100.     SVD( H, Q, D, R, TRUE, TRUE );
  101.     Q^.Show( 'S' );
  102.     D^.Show( 'D' );
  103.     R^.Show( 'R' );
  104.     S := matequals( S, mult( mult( Q, Diag( D ) ), Tran( R ) ) );
  105.     S^.Show( 'SDR`' );
  106.     
  107.     S := matequals( S, Ginv( H ) );
  108.     S^.Show( 'Ginv(H)' );
  109.     S := matequals( S, mult( S, H ) );
  110.     S^.Show( 'h*Ginv(h)' );
  111.     
  112.     dispose( h, killvmatrix );
  113.     dispose( S, killvmatrix );
  114.     dispose( Q, killvmatrix );
  115.     dispose( D, killvmatrix );
  116.     dispose( R, killvmatrix );
  117.     
  118.     Dispatch^.Declevel;
  119. End;
  120.  
  121. Procedure testShape;
  122. Var
  123.   H, S, Q, D, R : vmatrixptr;
  124. Begin
  125.     Dispatch^.Inclevel;
  126.     new( H, makematrix( 1, 1 ) );
  127.     new( S, makematrix( 1, 1 ) );
  128.     
  129.     h := matequals( h, Add( Ident( 5 ), Fill( 5, 5, 1 ) ) );
  130.     s := matequals( s, Vec( h ) );
  131.     S^.Show( 'Vec(h)' );
  132.     
  133.     s := matequals( s, shape( s, 5 ) );
  134.     S^.Show( 'Shape(Vec(s)),5' );
  135.     
  136.     s := matequals( s, Diag( h ) );
  137.     s^.show( 'Diagonal elements of h' );
  138.     Dispatch^.Declevel;
  139. End;
  140.  
  141. Procedure testFFT;
  142. Var
  143.   i,n : integer;
  144.   omega, sigma, xbar : double;
  145.   test,f,t, average, centered, power_spectrum : vmatrixptr;
  146.   covar, periodogram : vmatrixptr;
  147. Begin
  148.     n := 30;                           { n=2*3*5 }
  149.     omega := pi / n;
  150.     
  151.     new( t, makematrix( n, 1 ) );
  152.     new( f, makematrix( n, 2 ) );
  153.     new( test, makematrix( n, 2 ) );
  154.     new( average, makematrix( 1, 1 ) );
  155.     new( centered, makematrix( 1, 1 ) );
  156.     new( power_spectrum, makematrix( 2 * n, 1 ) );
  157.     new( covar, makematrix( 2 * n, 1 ) );
  158.     new( periodogram, makematrix( n, 1 ) );
  159.     
  160.     For i := 0 To n - 1 Do
  161.         t^.mm( i + 1, 1 )^ := sin( i * omega ) + cos( 0.3 * i * omega );
  162.     
  163.     f := matequals( f, fft( t, TRUE ) );
  164.     f^.show( 'time to frequency transform' );
  165.     
  166.     test := matequals( test, ch( t, fft( f, FALSE ) ) );
  167.     test^.show( 'frequency to time transform' );
  168.     
  169.     { get power spectrum, and serial correlations }
  170.     average := matequals( average, mult( Fill( 1, t^.r, 1 ), t ) );
  171.     xbar := average^.m( 1, 1 ) / n;
  172.     centered := matequals( centered, sub( t, Fill( t^.r, 1, xbar ) ) );
  173.     centered := matequals( centered, cv( t, Fill( t^.r, 1, 0 ) ) );
  174.     f := matequals( f, fft( centered, TRUE ) );
  175.     
  176.     For i := 1 To 2 * n Do
  177.         power_spectrum^.mm( i, 1 )^ := (sqr( f^.m( i, 1 ) ) + sqr( f^.m( i, 2 ) ) ) / n;
  178.     power_spectrum^.show( 'power spectrum' );
  179.     
  180.     covar := matequals( covar, fft( power_spectrum, FALSE ) );
  181.     sigma := covar^.m( 1, 1 );
  182.     covar := matequals( covar, submat( covar, 1, n, 1, 1 ) );
  183.     For i := 1 To n Do Begin  
  184.         covar^.mm( i, 1 )^ := covar^.m( i, 1 ) / sigma;
  185.         periodogram^.mm( i, 1 )^ := power_spectrum^.m( i, 1 ) / sigma;
  186.     End;
  187.     covar^.show( 'serial correlations' );
  188.     periodogram := matequals( periodogram, Mlog( periodogram ) );
  189.     periodogram^.show( 'periodogram' );
  190.     
  191.     dispose( t, killvmatrix );
  192.     dispose( f, killvmatrix );
  193.     dispose( test, killvmatrix );
  194.     dispose( average, killvmatrix );
  195.     dispose( centered, killvmatrix );
  196.     dispose( power_spectrum, killvmatrix );
  197.     dispose( covar, killvmatrix );
  198.     dispose( periodogram, killvmatrix );
  199.     
  200. End;
  201.  
  202. Procedure TestSums;
  203. Var
  204.   a, b : vmatrixptr;
  205. Begin
  206.     new( a, makematrix( 1, 1 ) );
  207.     new( b, makematrix( 1, 1 ) );
  208.     
  209.     a := matequals( a, Fill( 5, 5, 1 ) );
  210.     b := matequals( b, Sum( a, ALL ) );
  211.     b^.show( 'sum over all' );
  212.     b := matequals( b, Sum( a, ROWS ) );
  213.     b^.show( 'sum over rows' );
  214.     b := matequals( b, Sum( a, COLUMNS ) );
  215.     b^.show( 'sum over columns' );
  216.     
  217.     b := matequals( b, Sumsq( a, ALL ) );
  218.     b^.show( 'sum sq over all' );
  219.     b := matequals( b, Sumsq( a, ROWS ) );
  220.     b^.show( 'sum sq over rows' );
  221.     b := matequals( b, Sumsq( a, COLUMNS ) );
  222.     b^.show( 'sum sq over columns' );
  223.     
  224.     b := matequals( b, Cusum( a ) );
  225.     b^.show( 'Cusum of a' );
  226.     
  227.     a := matequals( a, b );
  228.     b := matequals( b, Mmax( a, ALL ) );
  229.     b^.show( 'max over all' );
  230.     b := matequals( b, Mmax( a, ROWS ) );
  231.     b^.show( 'max over rows' );
  232.     b := matequals( b, Mmax( a, COLUMNS ) );
  233.     b^.show( 'max over columns' );
  234.     
  235.     b := matequals( b, Mmin( a, ALL ) );
  236.     b^.show( 'min over all' );
  237.     b := matequals( b, Mmin( a, ROWS ) );
  238.     b^.show( 'min over rows' );
  239.     b := matequals( b, Mmin( a, COLUMNS ) );
  240.     b^.show( 'min over columns' );
  241.     
  242.     dispose( a, killvmatrix );
  243.     dispose( b, killvmatrix );
  244.     
  245. End;
  246. Procedure testElementary;
  247. Var
  248.   a : vmatrixptr;
  249. Begin
  250.     new( a, makematrix( 1, 1 ) );
  251.     a := matequals( a, Add( Ident( 5 ), Fill( 5, 5, 1 ) ) );
  252.     
  253.     Crow( a, 2, 3.0 );
  254.     a^.show( 'r2*3.0' );
  255.     Srow( a, 1, 2 );
  256.     a^.show( 'r1 <-> r2' );
  257.     Lrow( a, 1, 2, - 3.0 );
  258.     a^.show( '-3.0*r2 + r1' );
  259.     
  260.     Ccol( a, 2, 3.0 );
  261.     a^.show( 'c2*3.0' );
  262.     Scol( a, 1, 2 );
  263.     a^.show( 'c1 <-> c2' );
  264.     Lcol( a, 1, 2, - 3.0 );
  265.     a^.show( '-3.0*c2 + c1' );
  266.     
  267.     dispose( a, killvmatrix );
  268. End;
  269.  
  270. Begin
  271.     {main function}
  272.     testMfunctions;
  273.     testPatterned;
  274.     testDecomp;
  275.     testShape;
  276.     testFFT;
  277.     testSums;
  278.     testElementary;
  279. End.
  280.